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FAAH-like anandamide transporter (FLAT) regulates anandamide transport for hydrolysis and may be an 
attractive drug target for pain regulation. We aimed to discover potential FLAT antagonists from traditional 
Chinese medicine (TCM) using virtual screening, ligand-based drug design and molecular dynamics 
simulation (MD). Guineensine and Retrofractamide A exhibited high Dock Scores in FLAT. Consensus 
from multiple linear regression (MLR; R^ = 08973) and support vector machine (SVM; R^ = 0.7988) showed 
similar bioactivities for Guineensine and the FAAH-1 inhibitor (9Z)-l-(5-pyridin-2-yl-l,3>4-oxadiazol- 
2-yl)octadec-9-en-l-one. Contour of Guineensine to CoMFA and CoMSIA features also imply bioactivity. 
MD revealed shake or vibration in the secondary structure of FLAT complexed with Guineensine and 
(9Z)-l-(5-pyridin-2-yl-l,3)4-oxadiazol-2-yl)octadec-9-en-l-one. Ligand movement might contribute to 
protein changes leading to vibration patterns. Violent vibrations leading to an overall decrease in FLAT 
function could be the underlying mechanism for Guineensine. Here we suggest Guineensine as a drug-like 
compound with potential application in relieving neuropathic pain by inhibiting FLAT. 

Neuropathic pain is a multifactor neurogenic disorder caused by physical damage of neurons, cancer and 
other diseases \ Patients of neuropathic pain suffer from chronic pain as well as mental illnesses such as 
depression, anxiety and sleeping disorders^. To date, the mechanism of neuropathic pain remains 
unclear, making diagnosis and treatment difficult^'^'^. 

Anandamide is an endogenous cannabinoid formed by the N-acyl-phosphatidylethanolamine-selective phos- 
phodiesterase (NAPE-PLD) catalyzed hydrolysis of N-arachidonoyl-phosphatidyl-ethanolamine (NAPE)^, and 
has important physiological roles in pain regulation^. However, activity period of anandamide is short due to the 
rapid inactivation of anandamide by fatty acid amide hydrolyase (FAAH-1 )^'^. Catoblism of anandamide is 
associated with many different diseases, including cancer, cardiovascular disease, obesity, and particularly neuro- 
pathic pain^"^^. One emerging approach in controlling pain is the modulation of anandamide degradation by 
targeting FAAH-1 Several antagonists of FAAH have been successfully developed^^"^^. Recent findings 
suggest FAAH-1 cytosolic variant FAAH-like anandamide transporter (FLAT)^^ as a possible target for regulating 
pain. Decreased transportation of anandamide to FAAH-1 by inhibiting FLAT may be an alternative to direct 
antagonism of FAAH. 

In this study, we screen for drug-like compound against FLAT from TCM Database@Taiwan^°. Ligand based 
drug design methods were employed to predict bioactivity of the selected ligands. Molecular dynamics were 
employed to investigate underlying molecular mechanisms that may contribute to FLAT inhibition. 

Results 

Homology modeling and molecular docking. Suitability of rat proteins as templates for modeling human 
proteins was assessed by sequence alignment. Alignment of native rat FAAH-1 and native human FAAH-1 
sequences showed 79.7% identity and 89.8% similarity. Re-alignment following removal of ot2-helices (T9-T76) 
(termed FLAT sequences for clarification purposes) increased sequence identity and similarity to 86.1% and 
95.6%, respectively (Figure 3). We proceeded to model human FLAT structure using rat FLAT structure based on 
the high sequence identity and similarity of the FLAT sequences. Structural correctness of the modeled human 
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Figure 1 | Experimental procedure and structural basis of FLAT 
simulation, (a) Simplified scheme of experimental procedures, 
(b) Structural basis for FLAT structure simulation using FAAH-1. The a2- 
interacting loop (K255-L278; red) is the binding site opening loop, and the 
helices (P411-N435) colored in cyan are regions in FAAH-1 that interact 
with the membrane. Presence of the a2-helix (T9-T76; orange) in FAAH-1 
was the primary structural difference from FLAT. Human FLAT was 
modeled from rat FLAT structure, which was computationally prepared by 
deleting the a2 -helix region (amino acids T9-T76) in rat FAAH-1. 
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Figure 2 | Cartoon representation of the anandamide binding site and 
docking poses of TCM candidates within the binding site, (a) Enlarged 
view of the docking site (green) within the modeled human FLAT protein. 
The front and back sides of the binding site are depicted in red and purple, 
respectively, (b) Front view of docking site with docked ligands. (9Z)-l-(5- 
pyridin-2-yl- 1 ,3,4-oxadiazol-2-yl)octadec-9-en- 1 -one (control) , Guineensine, 
and Retrofractamide A are shown in orange, blue, and green, respectively. 
Ser217 and Ile238 are amino acids found adjacent to the binding site, (c) Side 
view of docking site with docked ligands (45 degrees relative to b). 
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Figure 3 | Sequence alignment of target human FLAT sequence with template rat FLAT sequence. FLAT sequences refer to that of rat FAAH-1 (PDB: 
3K84) and human FAAH-1 (SwissProt: 000519) in which amino acids T9-T76 have been removed. Sequence identity and similarity were 86.1% and 
95.6%, respectively. 
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Figure 4 | Ramachandran validation of the modeled human FLAT structure. A total of 98.2% of the residues were distributed in the favored region. The 
remaining 1.8% were located in the allowed region. No residues have psi or phi angles in the disfavored regions. 



FLAT structure was evaluated with the Ramachandran plot. A total 
of 491 residues (98.2%) were distributed in the favored region 
(Figure 4). Table 1 lists the nine residues (1.8%) distributed in the 
allowed region. Results of the Ramachandran plot suggest that the 
modeled human FLAT structure is correct. 



Structure scaffolds of the control and TCM candidates are shown 
in Figure 5. Candidate selection was primarily based on structural 
similarities to the control and Dock Score which considers the in- 
ternal energies of ligands and their corresponding protein-ligand 
interaction energy. Guineensine^^ and Retrofractamide A^^, 
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Table 1 | Residues of the modeled human FLAT protein structure 
with (p and v|/ angles located within the allowed region 



Residue 


Amino acid 


Region 


Coordinate 


90 


Arg 


Allowed 


( 80.28, 9.05) 


136 


Tyr 


Allowed 


(-38.98, 114.04) 


153 


Ser 


Allowed 


(-125.42,81.78) 


191 


Ser 


Allowed 


(73.17, -4.14) 


240 


Gly 


Allowed 


(-148.28, 11.08) 


312 


Pro 


Allowed 


(-70.50, 95.67) 


335 


Asn 


Allowed 


(74.14, -2.62) 


484 


Ala 


Allowed 


(-79.76, -38.61) 


488 


Ala 


Allowed 


(-112.69, 80.14) 



compounds originating from Piper longum^^'^^ have Dock Scores of 
61.45 and 54.01, respectively, and were selected as candidates 
(Table 2). Front and side views of the docking poses of the TCM 
candidates and the control at the modeled anandamide binding site 
are illustrated in Figure 2b and 2c. Ligand-protein interactions can be 
extrapolated from Figure 6. For the control, 16 hydrophobic contacts 
were formed with 14 residues in FLAT (Figure 6a). Hydrogen bonds 
(H-bonds) were formed with Ser218, Ile239, Gly240, and Ser242. 
Guineensine was stabilized by solely hydrophobic interactions 
(Figure 6b). Similar to the control, Retrofractamide A formed both 
H-bonds and hydrophobic interactions with FLAT. Since Guineensine 
and Retrofractamide A primarily differ in the number of carbons in 
the hydrophobic tail (Figure 5b), the difference in hydrophobic tail 
length may contribute to differences in scoring functions (Table 2). 
Based on Ligplot analysis, the O atom located on the hydrocarbon 
chain was situated to form H-bonds with Metl92, Serl94, Ser242 in 
Retrofractamide A. By contrast, the longer hydrocarbon chain in 
Guineensine pushes the same O atom away from the H-bond forming 
amino acids, resulting in Guineensine being stabilized by only hydro- 
phobic interactions (Figure 6b). Cross comparing interaction residues 
among different test ligands reveals that Metl92, Phel93, Serl94, 
Ile239, and Val492 are residues with which all test ligands formed 
interactions with (Table 3). The residues functioned to anchor the 
hydrocarbon tails of the ligands (Figure 7). 




Retrofractamide A 
Guineensine 




Figure 5 | Scaffold structures of control compound and TCM candidates. 

(a) the alpha-ketoheterocycle inhibitor (9Z)-l-(5-pyridin-2-yl-l,3,4- 
oxadiazol-2-yl)octadec-9-en-l-one (control), and TCM candidates (b) 
Guineensine and Retrofractamide A. Guineensine and Retrofractamide A 
differ in the number of carbons found in the hydrophobic tail region. The 
number of carbons for Guineensine is depicted in blue and that for 
Retrofractamide A is depicted in gray. 
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Figure 6 | Protein-ligand interactions analyzed by Ligplot. (a) (9Z)-l-(5-pyridin-2-yl-l,3>4-oxadiazol-2-yl)octadec-9-en-l-one, (b) Guineensine, and 
(c) Retrofractamide A. 



Activity prediction using MLR and SVM. Representative descrip- 
tors determined by GFA and used for constructing MLR and SVM 
models were: C_Count, ES_Sum_dsCH, ES_Sum_aasQ ES_Sum_ 
ssNH, ES_Sum_dsN, ES_Sum_aaN, ES_Sum_dO, ES_Sum_sF, 
AverageBondLength, Num_Chains. Descriptions of each property 
are listed in Supplementary Table SL In general, antagonistic 
activity against FLAT was associated with carbon atom numbers, 
electrotopolocial state, atom bond length, and unbranched chain 
numbers. The optimum GFA model with an R^ of 0.9403 was: 
GFATempModelA = 62.509 + 59.957 x C_Count - 210.91 x £S 

_Sum_dsCH + 37. 751 xES_Sum_aasC- 432.98 xES_Sum 

-SsNH + 60.838 xES-Sum-dsN + 5.0225 xES-SumMaN 

+ 87 .622 X ES-Sum-dO-2.4797 X ES-Sum^F -258.78 x 

AverageBondLength — 30.318 x Num_Chains 

The coefficient of determination (R^) for the MLR and SVM mod- 
els constructed using the aforementioned descriptors are 0.8973 
(Figure 8a) and 0.7988 (Figure 8b), respectively. Correlation between 
observed and predicted bioactivity values indicate good prediction 
accuracy of both models. Predicted pICso of the control and TCM 
candidates using the MLR and SVM models are listed in Table 2. 
Both models predict bioactivity for the TCM candidates. Similar 
bioactivities ranging between 7.72-7.89 were predicted for the three 
compounds with MLR. SVM model predictions show Guineen- 
sine and (9Z)-l-(5-pyridin-2-yl-l,3,4-oxadiazol-2-yl)octadec-9-en- 
1-one to have bioactivities approximately 1.5 logs higher than 
Retrofractamide A. To reduce possible bias caused by one single 
model, we used a simplified consensus to summarize our bioactivity 
prediction results. The criteria for each score given is described in 
Table 2. The combined results of our MLR and SVM models indicate 
that Guineensine and the control have similar predicted bioactivities 
while Retrofractamide A may have lower bioactivities (Table 2). The 
variation in bioactivity could be related to the number of carbons 
within the hydrophobic tail since carbon numbers were the primary 
difference between Guineensine and Retrofractamide A. GFA results 



in which the number of C atoms is a representative descriptor of 
bioactivity support this view. 

Activity prediction using CoMFA and CoMSIA. PLS results of 
CoMFA and CoMSIA models are tabulated in Table 4. Steric field 
was the primary determinant of bioactivity in the CoMFA model. At 
an optimal number of components (ONC) of 6, the cross validation 
correlation coefficient (q^) of 0.690 and non-cross-validation 
correlation coefficient (r^) of 0.954 indicated a confident model. 
CoMSIA models considered multiple factors, and the possible 
models at an ONC of 6 are shown in Table 4. The CoMSIA model 
integrating steric, hydrophobic, and H-bond acceptor properties 
exhibited the highest cross validation (q^ = 0.802) and non-cross 
vaUdation (r^ = 0.948) values, thus was selected as the CoMSIA 
model of choice. 

Correlation between the observed and predicted bioactivity values 
generated by the CoMFA and CoMSIA are shown in Figure 8c and 
8d, respectively. Correlation coefficients (R^) of 0.9542 for the 
CoMFA model and 0.9478 for the CoMSIA model indicated good 
prediction strength of both models. The range of residuals between 
predicted and observed value is < ±1.0 (Table 5). Superimposition 
of aligned ligands with CoMFA feature contours indicate that ligand 
aromatic heads and bends in the hydrophobic tail region contoured 
to the steric favoring bulk (yellow) in the CoMFA map (Figure 9a). 
Additionally, the hydrophobic tail matches the contour disfavoring 
steric bulk (green). Alignment of the ligands with the bioactivity 
feature maps suggest activity of the ligands against the target protein 
FLAT. As expected, the CoMSIA feature map (Figure 9b) differed 
from CoMFA. Steric favoring contours (yellow), as well as steric 
disfavoring contours (green), were located near the hydrophobic tail. 
Both hydrophilic (cyan) and hydrophobic (gray) contours were 
located near the aromatic head. From the contour of the aromatic 
head region to the hydrophobic favoring bulk, we speculate that the 
ligands should exhibit some bioactivity despite the lack of bulk at the 
steric favoring region. Results from CoMFA and CoMSIA are con- 
sistent with results from MLR and SVM predictions, and consistently 
suggest bioactivity of the TCM candidates toward FLAT. 



Table 3 | Analysis of residues in FLAT that interact with test ligands 

AminoAcid MFSS IGSGCVYQLYLE LFVMW 



Sequence 192 193 194 218 239 240 242 269 270 271 272 274 279 336 373 374 381 382 492 496 532 
Control OOOOOOO-OO - - - OOO-OOOO 

Guineensine 000-0-----0-----0-0- 
Retrofractamide A 000-0-00---00-----0- 



O: Residue which interact with ligand 
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Figure 7 | Spatial orientation of key interaction residues with regard to test ligands. (a) (9Z)-l-(5-pyridin-2-yl-l,3>4-oxadiazol-2-yl)octadec-9-en-l- 
one (mint), (b) Guineensine (blue), and (c) Retrofractamide A (violet). Common residues are illustrated by yellow surfaces. 



Molecular dynamics (MD) simulation. Molecular dynamics 
simulations were employed to investigate dynamic stability of the 
complexes. Characteristics of the protein-ligand complexes are 
summarized in Table 6. Complexes of the TCM candidates 
Guineensine and Retrofractamide A had lower total energy levels 
than the control complex. Globally speaking, binding of TCM 
candidates did not significantly alter protein characteristics 



compared to the control FAAH-1 antagonist. On a local scale, 
calculated RMSD values and total drifts (Table 6) indicate that 
Retrofractamide A complex was most stable, followed by (9Z)-1- 
(5-pyridin-2-yl- 1 ,3,4-oxadiazol-2-yl)octadec-9-en- 1 -one. Guineen- 
sine was the least stable of the three complexes. 

DSSP analysis allowed visualization of structural changes that 
contributed to stability differences. Three distinctive patterns, (1) 
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Figure 8 | Correlation plot of observed bioactivity verses predicted bioactivities generated by different models, (a) MLR (b) SVM (c) CoMFA (d) 
CoMSIA. 
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U.U4 


n 1 o 

U. 1 V 


SADH 


6 


0.73 


0.95 


0.28 


101.24 


0.35 




0.48 


0.02 


0.14 


SADE 


6 


0.71 


0.93 


0.32 


78.44 


0.73 


0.00 




0.07 


0.21 


SAEH 


6 


0.74 


0.95 


0.28 


106.00 


0.36 


0.00 


0.49 




0.14 


SDEH 


6 


0.69 


0.94 


0.29 


94.80 


0.42 


0.00 


0.56 


0.01 




ADEH 


6 


0.76 


0.94 


0.94 


93.90 




0.00 


0.78 


0.04 


0.19 


SADEH 


6 


0.80 


0.95 


0.28 


101.24 


0.35 


0.00 


0.48 


0.02 


0.14 



ONC: Optimal number of components 
S: Steric 

SEE: Standard error of estimate 

E: Electrostatic 

F: F-test value 

H: Hydrophobic 

PLS: partial least squares 

D: Donor 

*: Optimum prediction model 
A: Acceptor 



without shake or vibration, (2) shake, or (3) vibration, were observed 
in the three complexes (Figure 10). The non-vibrating pattern 
impUes a stable secondary protein structure of small changes. The 
shake pattern refers to isolated peaks of changes whereas the vibra- 
tion pattern exhibits constant vibrations over a period of time. The 
shake pattern can be viewed as a pre-vibration signal. As shown in 
Figure 10a, all three patterns are present in the control complex. 
Isolated shaking patterns were present from 15-26 ns, and accumu- 
lated to continuous vibrations from 26-40 ns. Guineensine complex 
was stable up to 6 ns, but exhibited vibration patterns after 6 ns to 
the end of MD (Figure 10b). Retrofractamide A complex remained 
stable during MD, with minor shakes observed in residues 497-577 
after 30 ns (Figure 10c). The fluctuations observed here may provide 
an explanation for the larger RMSDs and drifts recorded in Table 6. 
Root mean square fluctuation (RMSF) analysis on docking inter- 
action residues provide information on residues that may be import- 
ant for the protein DSSP changes and ultimately bioactivity against 
FLAT. Large fluctuations at Ile239, Tyr272, and Gln274 could induce 
ligand instability, and in turn affect the stability of FLAT. (Figure 11, 
Table 7). 



Several analyses were conducted to determine the effects of 
TCM ligands on FLAT structure with (9Z)-l-(5-pyridin-2-yl-l, 
3,4-oxadiazol-2-yl)octadec-9-en-l-one as the point of reference 
(Figure 12). Solvent accessible surface area (SASA) remained rela- 
tively constant for protein complexes (Figure 12a) and individual 
ligands (Figure 12b), indicating no significant changes in protein 
folding. Average ligand SASAs of control, Guineensine and Retro- 
fractamide A are 6.895 nm/NS^ 6.458 nm/NS" and 5.488 nm/NS^ 
respectively, which is most likely related to the hydrophobic tail 
lengths of each ligand. Changes in complex compactness were eval- 
uated based on radius of gyration (R^) which measures the mass of 
atoms relative to the center of mass of the complex. Average and 
maximum levels of R^ for each complex are tabulate in Table 8. Over 
the course of time, an increase in R^compiex starting at approximately 
28 ns was recorded for the control, and maximum R^compiex of 
4.184 nm was reached at 38.46 ns. In Guineensine, changes from 
the baseline trajectory were observed around 10 ns and peaked to 
R^compiex of 4.7 nm at 18.32 ns. The R^compiex for Retrofractamide A 
was relatively static. In agreement with R^compiex results, higher 
%iigand, was recorded for Guineensine compared to the other two 
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Table 5 | Bioactivity residuals^ between observed and predicted bioactivity values generated with the CoMFA and CoMSIA models 

CoMFA CoMSIA 



Compound Index^ Observed Bioactivity^ (PIC50) 


Predicted PIC50 


Residual 


Predicted pl5o 


Residua 


4 


6.47 


6.56 


-0.09 


6.50 


-0.03 


5 


7.55 


7.53 


0.02 


7.60 


-0.05 


6 


7.72 


7.73 


-0.01 


7.79 


-0.07 


8 


7.82 


7.85 


-0.03 


7.87 


-0.05 


9 


7.89 


7.91 


-0.02 


7.86 


0.03 


22 


4.96 


5.00 


-0.04 


5.14 


-0.1 8 


23 


5.20 


4.81 


0.39 


4.73 


0.47 


25 


6.14 


6.1 9 


-0.05 


6.05 


0.09 


26 


6.68 


6.71 


-0.03 


6.68 


0 


27 


7.55 


7.39 


0.1 6 


7.42 


0.1 3 


29 


7.74 


7.90 


-0.1 6 


7.87 


-0.1 3 


31 


8.62 


8.63 


-0.01 


8.46 


0.1 6 


32 


8.68 


8.79 


-0.1 1 


8.59 


0.09 


33 


8.66 


8.53 


0.1 3 


8.41 


0.25 


35 


8.30 


8.46 


-0.16 


8.42 


-0.1 2 


36 


6.54 


6.88 


-0.34 


6.90 


-0.36 


37 


7.09 


7.28 


-0.1 9 


7.20 


-0.1 1 


38 


8.03 


7.69 


0.34 


7.76 


0.27 


41 


8.10 


8.18 


-0.08 


8.22 


-0.1 2 


42 


8.07 


8.43 


-0.36 


8.65 


-0.58 


44 


9.60 


9.38 


0.22 


9.25 


0.35 


45 


9.46 


9.46 


0 


9.43 


0.03 


46 


9.35 


9.50 


-0.15 


9.34 


0.01 


47 


8.62 


8.87 


-0.25 


8.84 


-0.22 


48 


7.34 


7.32 


0.02 


7.43 


-0.09 


49 


7.89 


7.69 


0.20 


7.97 


-0.08 


50 


8.66 


7.93 


0.73 


8.02 


0.64 


51 


8.74 


8.50 


0.24 


8.46 


0.28 


52 


8.19 


8.19 


0 


8.41 


-0.22 


53 


8.64 


8.54 


0.1 0 


8.66 


-0.02 


57 


9.40 


9.44 


-0.04 


9.32 


0.08 


58 


8.68 


8.91 


-0.23 


9.06 


-0.38 


59 


6.82 


7.31 


-0.49 


7.1 5 


-0.33 


60 


7.54 


7.76 


-0.22 


7.69 


-0.1 5 


61 


7.85 


7.94 


-0.09 


7.73 


0.1 2 


62 


8.21 


8.37 


-0.1 6 


8.49 


-0.28 


63 


8.08 


8.22 


-0.14 


8.1 8 


-0.1 0 


64 


8.64 


8.44 


0.20 


8.77 


-0.1 3 


65 


9.10 


9.09 


0.01 


9.07 


0.03 


66 


9.19 


9.37 


-0.1 8 


9.36 


-0.1 7 


68 


9.82 


9.50 


0.32 


9.47 


0.35 


69 


9.46 


8.90 


0.56 


8.84 


0.62 


28* 


7.55 


7.07 


0.48 


6.85 


0.70 


30* 


8.08 


7.96 


0.12 


7.78 


0.30 


34* 


8.31 


8.75 


-0.44 


8.57 


-0.26 


39* 


8.07 


7.91 


0.16 


7.84 


0.23 


40* 


8.46 


8.28 


0.18 


8.35 


0.1 1 


43* 


9.37 


9.09 


0.28 


9.13 


0.24 


54* 


9.57 


9.12 


0.45 


9.35 


0.22 


55* 


9.70 


9.27 


0.43 


9.30 


0.40 


56* 


9.60 


9.48 


0.12 


9.63 


-0.03 


67* 


9.46 


9.50 


-0.04 


9.40 


0.06 


^Residual = observed pICso-predicted pICso- 

^Ligands and corresponding experimental bioactivity values v/ere adopted from [30; Sit 201 0]. 



*: test set 



test ligands (Figure 12d). Intriguingly, R^compiex trajectories seemed 
to be related to pattern changes in Figure 10. Changes in R^compiex 
coincided with shakes or vibrations observed for the control and 
Guineensine. 

The time dependent movement of atoms from their initial posi- 
tions (MSD) were in agreement with gyration results. An increase in 
MSDcompiex (indicated by the increase in slope) was observed for the 
control (Figure 12e) as the complex becomes more relaxed after 
25 ns (Figure 12c). A two-stage change in MSDcompiex was observed 



for Guineensine (Figure 12e). The steep slope in stage 1 (0-17 ns) 
matched the expansion of Guineensine complex (Figure 12c) and 
structure vibration (Figure 10b). During stage 2, the decrease in 
MSDcompiex slope Coincided with increased Guineensine compact- 
ness and lower intensity vibrations. Displacement of ligands may also 
contribute to protein compactness and vibrations. Trajectory of con- 
trol MSDiigand was similar to that for Rgcomplex. For Guineensine, 
the MSDiigand was similar to the trend of its MSDpj-otein- It is possible 
that in both cases, ligand shift was a primary force in initiating 
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Steric (disfavored) 



hydrophobic yWi ' 




iSteric (disfavored) 




Steric (favored) 



Figure 9 | Contour of superimposed ligands to CoMFA and CoMSIA 
feature maps, (a) The aromatic head and bend in the hydrophobic tail 
contoured to the steric favoring region (yellow), and the hydrophobic tail 
regions contoured to the steric disfavoring region in the CoMFA map. 
(b) In the CoMSIA map, ligand tail regions contoured to both the steric 
favoring (yellow) and disfavoring (green) regions. Aromatic heads of the 
ligands contoured to both the hydrophobic (gray) and hydrophilic (cyan) 
features. (9Z) - 1 - (5-pyridin-2-yl- 1 ,3,4-oxadiazol-2-yl) octadec-9-en- 1 - 
one, Guineensine, and Retrofractamide A (green) are represented in blue, 
red, and green, respectively. 

protein changes recorded in Figure 10 and Figure 12c. The high 
RMSDpj-otein of the control as opposed to the TCM candidates 
(Figure 12g) suggest that complexes formed by the TCM candidates 
cause less structural perturbation to the FLAT protein. No significant 
differences in RMSDngand was observed among the three compound 
(Figure 12h). Distance matrices of the three candidates v^ere also 
similar (Figure 13). 

Discussion 

The selected TCM compounds Guineensine and Retrofractamide A 
are near identical chemical structures with the only difference being 
the length of the hydrocarbon chain length (Figure 5). The longer tail 
of Guineensine enables its stabilization by hydrophobic interactions 
within the FLAT protein (Figure 6). By contrast, the shorter tail in 
Retrofractamide A formed H-bonds with Metl92, Serl94, Ser242 
within the hydrophobic interior of the protein. Residues Met 192, 
Phel83, Serl94, Ile239, Val492 were important contact residues that 
stabilized the tail regions within the docking site. Bioactivity results 
obtained through QSAR model predictions also suggested an asso- 
ciation between bioactivity strength and chain length. Guineensine 
and the control were both predicted to have higher bioactivities than 
the shorter Retrofractamide A. MD results further revealed protein 
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Without Vibration 




20 
Time (ns; 



20 

Time (ns) 



□ coil 



p-sheet 



p-bridge 



bend 



□ turn 



a-helix 



3-helix 



Figure 10 | Secondary structure changes observed during the 40 ns MD simulation, (a) (9Z)-l-(5-pyridin-2-yl-l,3>4-oxadiazol-2-yl)octadec-9-en-l- 
one, (b) Guineensine, and (c) Retrofractamide A. The three complexes demonstrate distinctive change patterns. The "shake" pattern refers to a sudden 
peak of vibration; the "vibration" pattern refers to observation of continuous vibrations. 



W532 
M496 
V492 
F382 
L381 
E374 
L373 
Y336 
L279 
|Q274 
w Y272 
q: V271 
0270 
G269 
S242 
G240 
1239 
S218 
S194 
F193 
M192 





T 



I I I I I I I I I I I I I I I I I I I 



0.00 0.02 0.04 0.06 0.08 0.10 
RMSF(nm) 

Figure 1 1 | Root mean square analysis of FLAT residues implicated in 
docking. Yellow, red, and blue bars represent data for (9Z)-l-(5-pyridin- 
2-yl-l,3,4-oxadiazol-2-yl)octadec-9-en-l-one, (b) Guineensine, and (c) 
Retrofractamide A, respectively. 



secondary structure changes in Guineensine and the control that 
were not observed in Retrofractamide A. Since ability to initiate 
unstable vibrations in FLAT showed similar trends to bioactivity 
(Figure 10), we propose that a possible link might exist between these 
characteristics. Large fluctuations at Ile239, Tyr272, and Gln274 
could contribute to bioactivity by inducing ligand instability, chan- 
ging protein compactness and affecting the stability of FLAT. Based 
on the combined results from docking, QSAR predictions, and MD 
simulation, Guineensine exhibits similar patterns to the control and 
could be a likely candidate for further studies. 



Table 7 


Root mean square 


fluctuations (RMSF) 


(nm) recorded for 


residues involved in docking 










RMSD (nm) 




Residue 


Control 


Guineensine 


Retrofractamide A 


Ml 92 


0.0435 


0.054 


0.0671 


F193 


0.0232 


0.038 


0.0784 


5194 


0.0136 


0.0123 


0.0164 


5218 


0.0098 


0.0242 


0.0433 


1239 


0.0723 


0.0716 


0.0468 


G240 


0.012 


0.0123 


0.01 12 


5242 


0.01 12 


0.0253 


0.0142 


G269 


0.012 


0.0363 


0.0098 


C270 


0.0385 


0.0552 


0.0387 


V271 


0.0286 


0.0693 


0.0767 


Y272 


0.0825 


0.0745 


0.0256 


Q274 


0.0481 


0.0812 


0.0235 


L279 


0.0664 


0.0636 


0.0757 


Y336 


0.0738 


0.0197 


0.0208 


L373 


0.0278 


0.0266 


0.0331 


E374 


0.068 


0.055 


0.0604 


L381 


0.0589 


0.0806 


0.07 


F382 


0.0277 


0.0188 


0.02 


V492 


0.0462 


0.0283 


0.028 


M496 


0.0477 


0.0647 


0.052 


W532 


0.0187 


0.0235 


0.0205 
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Figure 12 | Analysis of MD trajectories generated by Gromacs. Trajectories for (a) protein solvent accessible surface area (SASA), (b) ligand SASA, 
(c) protein radius of gyration (Rg), (d) ligand R^, (e) protein mean square deviation (MSD), (f) ligand MSD, (g) protein root mean square deviation 
(RMSD), and (h) ligand RMSD are shown. Data for (9Z)-l-(5-pyridin-2-yl-l,3)4-oxadiazol-2-yl)octadec-9-en-l-one, guineensine, and Retrofractamide 
A are shown in mint, blue, and violet, respectively. 
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Table 8 | Radius of gyration (Rg) calculated for the control and TCM candidates 

Radius of Gyration (Rg; nm) 



ControP Guineensine Retrofractamide A 





Average 


Maximum 


Average 


Maximum 


Average 


Maximum 


Protein 


2.596 


4.184 


3.100 


4.700 


2.170 


2.292 


Mainchain 


2.562 


4.153 


3.152 


4.794 


2.163 


2.299 


Sidechain 


2.602 


4.188 


3.093 


4.685 


2.172 


2.291 


Backbone 


2.556 


4.160 


3.148 


4.792 


2.157 


2.289 


C alpha 


2.561 


4.165 


3.153 


4.792 


2.162 


2.292 


Ligand 


0.623 


0.720 


1.621 


5.054 


0.569 


0.659 



^Control: (9Z)-1 -(5-pyridin-2-yl-l ,3,4-oxadiazol-2-yl)octadec-9-en-l-one 



The TCM compound Guineensine was predicted to have similar 
bioactivities with the demonstrated FAAH-1 inhibitor (9Z)-l-(5- 
pyridin-2-yl-l,3,4-oxadiazol-2-yl)octadec-9-en-l-one. Molecular in- 
sights suggest that the inhibitory mechanism may be related to the 
ability to initiate unstable vibrations in FLAT. Guineensine induced 
violent vibrations in the FLAT complex. (9Z)-l-(5-pyridin-2-yl- 
l,3,4-oxadiazol-2-yl)octadec-9-en-l-one induced similar, albeit less 
extreme, vibration. Gromacs analysis results suggest that ligand 
movement within the binding site might be associated with protein 
compactness changes which further translates into global protein 
structure vibrations. Based on these observations, we suggest that 
Guineensine may be a possible ligand for inhibiting FLAT trans- 
portation of anandamine. These findings may have important impli- 
cations in designing better alternative for managing neuropathic 
pain. 

Methods 

Figure la summarizes the overall procedure used in this study. 

Homology modeling. The crystal structure of rat fatty-acid amide hydrolase 1 (PDB: 
3K84)^^ was downloaded from Protein Data Bank. For clarification purposes, all 
sequence numberings referred to follow those of rat FAAH-1 (PDB: 3K84). The a2 
helix (from T9 to T76) of FAAH-1 was computationally removed to form the FLAT 
structure of rats (Figure lb). The homology model of human FLAT was constructed 
by using rat FLAT as the template structure and human FAAH-1 (SwissProt Index: 
000519) as the template sequence. RAMPAGE^^ was used to verify validity of the 
predicted structure. 

Molecular docking. The modeled human FLAT structure was applied for molecular 
docking. Over 30,000 compounds were downloaded from TCM Database@Taiwan^*' 



and docked into the anandamide binding site of the modeled human FLAT structure 
(Figure 2a). All compounds were prepared with CHARMm (Chemistry at HARvard 
Molecular Mechanics)^^ prior to docking to add missing hydrogen atoms. The alpha- 
ketoheterocycle inhibitor (9Z)-l-(5-pyridin-2-yl-l,3,4-oxadiazol-2-yl)octadec-9-en- 
1-one in rat FAAH-1 (PDB: 3K84) was used as the controP\ Control and TCM 
compounds were docked using the LigandFit module in Discovery Studio 2.5 (DS 
2.5)^"^. For each protein-ligand complex, five different conformations were generated 
based on default DS 2.5 settings. Each generated conformation was compared to the 
conformation of (9Z)-l-(5-pyridin-2-yl-l,3,4-oxadiazol-2-yl)octadec-9-en-l-one 
within the rat fatty-acid amide hydrolase 1 (PDB:3K84) crystal structure. The 
conformation most similar to the control in 3K84 and with the highest Dock Score 
was selected. 

Activity prediction using quantitative structure activity relationship (QSAR) 
Models. A total of 52 ligands^^ were selected to calculate molecular properties using 
the Calculate Molecular Properties protocol packaged in DS 2.5. Genetic function 
approximation (GFA) was used to identify ten representative descriptors associated 
with bioactivity. The calculated models were ranked by square correlation coefficient 
(R^), and the highest R^ model was used to predict bioactivities for the training set and 
test set. The linear Multiple Linear Regression (MLR) model was constructed with 
MATLAB^^, and the non-liner model utilizing Support Vector Machine (SVM)^^'^*^ 
was constructed with LibSVM^''. 

Comparative force field analysis (CoMFA) and comparative similarity indices 
analysis (CoMSlA) were applied to construct 3D-QSAR models. Ligand alignment 
was completed using the atom-fit module of SYBYL-X 1.1^°. In CoMFA, steric and 
electrostatic fields were calculated using Lennard- Jones potential (LJP) and 
Coulombic potential, respectively. In CoMSIA, descriptors including steric, hydrogen 
bond acceptor and donor, hydrophobic and electrostatic fields were calculated by 
Gaussian functions. Partial least squares (PLS)^^ was used for property analysis. Cross 
validation model was ranked by cross-validated coefficient (q^), whereas non-cross 
validation was ranked by conventional correlation coefficient (r^), standard error of 
estimate (SEE) and standard error. The model with the highest q^, r^ values, and the 
lowest SEE and standard error was selected as the optimum model, and used to 
predict bioactivities of the TCM candidates. 
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Molecular dynamics (MD) simulation. All ligands are prepared by using web server 
SwissParam (http://swissparam.ch/)^^ prior to MD simulation. MD was performed 
using GROMACS 4.0.7^^. The force field applied for simulation is detailed 
elsewhere^'*. Ligands and human FLAT were combined and immersed into a buffer 
(or solution) containing cubic box at a buffering distance of 1.2 nm between the 
complex and the edge of the box. Sodium and chloride ion were added to neutralize 
complex charge. Complex energies were subsequently minimized with the Steepest 
Descent method for 5000 steps. The last frame of each energy- minimized structure 
was used as the initial frame for MD simulation. Electrostatic interactions were 
calculated with the Particle-Mesh Ewald (PME) method^^ The cutoff for PME was 
1.0 nm. The time step was set at 2 fs and the number of steps set to 20000000, 
accumulating to a total MD simulation time of 40 ns. 

MD trajectories were analyzed with built-in Gromacs tools. The secondary struc- 
ture database (DSSP) was installed into Gromacs to analyze protein secondary 
structure changes^^. Program g_energy was used to analyze potential energy, kinetic 
energy, total energy, temperature, pressure, volume, density, pV and enthalpy 
changes. Program g_gyrate was used the measure the radius of gyration. Program 
g_sas was used to compute interaction surface areas between solvent molecules and 
complexes. Program ^_m(imaf was used to generate distance matrices which calculate 
the smallest distance between each residue pairs. 
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